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We present a detailed study of the exchange-correlation hole and exchange-correlation energy 
per particle in the Si crystal as calculated by the Variational Monte Carlo method and predicted 
by various density functional models. Nonlocal density averaging methods prove to be success- 
ful in correcting severe errors in the local density approximation (LDA) at low densities where 
the density changes dramatically over the correlation length of the LDA hole, but fail to provide 
systematic improvements at higher densities where the effects of density inhomogeneity are more 
subtle. Exchange and correlation considered separately show a sensitivity to the nonlocal semicon- 
ductor crystal environment, particularly within the Si bond, which is not predicted by the nonlocal 
approaches based on density averaging. The exchange hole is well described by a bonding orbital 
picture, while the correlation hole has a significant component due to the polarization of the nearby 
bonds, which partially screens out the anisotropy in the exchange hole. 
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I. INTRODUCTION 

Density functional theory (DFT) is the leading theo- 
retical tool for the study of the material properties of 
solids. It is based on the characterization of the ground- 
state energy of the inhomogeneous electron g func- 
tional of the density, whose optimization condition can 
be expressed in terms of the self-consistent single-particle 
Kohn-Sham equation!!! The success of the method lies 
in the ease of use and surprising effectiveness of the 
standard approximation, the local density approxima- 
tion (LDA). This models the key component of the func- 
tional, the exchange-correlation energy E xc [n], and the 
exchange-correlation potential used in the Kohn-Sham 
equation with the assumption that the inhomogeneous 
electron gas at any location in space behaves like the ho- 
mogeneous electron gas at a density equal to the local 
density. As a result, the formidably complex and nonlo- 
cal relationship between the energy and the density due 
to interparticle correlations in the inhomogeneous envi- 
ronment can be converted into a local functional depen- 
dent on input from the comparatively simpler homoge- 
neous electron gas. This approximation has been quite 
successful at obtaining useful and suprisingly accurate 
estimates for bulk ground-state properties of many solids 
and the basic ground-state structure of solids, surfaces 
and molecules. 

However, there remains a significant need for the de- 
velopment of more accurate density functionals in several 
areas. Applications in molecular and solid-state chem- 
istry require the calculation of total energies to at least 
an order of magnitude greater precision than the LDA 



can provide. Excited state properties, in particular the 
band-gap of semiconductors, are not typically obtain- 
able with quantitative accuracy. In addition, in mate- 
rials in which electron correlations are important, the 
LDA often fails to predict qualitative ground-state prop- 
erties. There have been numerous attempts in recent 
years to develop exchange-correlation functionals beyond 
the LDA_ iiicluding generalized gradient approximations 
(GGA's)QoO'El0 based on gradient-expansions about the 
LDA, model nonlocal potentials,™ seifiinteraction cor- 
rections,^ hybrids with Hartree-Fockliil and configura- 
tion interactional or other map¥Tbj^^heories,E3o and 
orbital-dependent functionals JlalJ^IIZHIa These have led 
to significant improvements in accuracy of DFT calcu- 
lations; but a functional of systematically quantitative 
accuracy for many quantum chemistry and solid-state 
applications has yet to be achieved. 

A key relation that has been particularly fruitful in 
the development of densifaz-functional theory is the adia- 
batic connection formulaEj which relates the exchange- 
correlation energy to a coupling-constant integrated 
exchange-correlation hole. The exchange-correlation 
hole, n xc , measures the change in density from its mean 
value at each point in a system, given the observation 
of an electron at one particular point, providing a simple 
visualization of the effects of electron correlation in an in- 
homogeneous material. By integrating this quantity over 
a family of systems characterized by interaction coupling 
constant A and kept at fixed density, the kinetic energy 
cost to create the hole at full coupling, A=l, is incorpo- 
rated into it. As a result, the Coulomb interaction energy 
of the integrated hole reflects the total energy associated 
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with creating the hole. 

The adiabatic connection formula, by making explicit 
the relation between the exchange-correlation energy and 
the exchange-correlation hole has been the impetus, both 
explicitly and implicitly, to many proposed corrections 
to the LDA. The average density approximation!! (ADA) 
and weighted density approximation^ (WDA) make ex- 
plicit use of this relation to construct nonperturbative 
nonlocal density functionals. WDA in particular has at- 
tracted attention as a potential tool with applications 
for quantum chemistty£jJ-,pptallic surfaces^ and bulk 
solid-state structure.BlBBlHil It is apparently limited by 



a lack of consistency in its behaviorErE3 and typically 
requires the use of various context-dependent models or 
extensions to achieve.-pptimal results for different appli- 
cations]^ The PW91B and PBEQ forms of the GGA ex- 
plicitly build in a number of exact and approximate prop- 
erties of the system average of the exchange-corrplation 
hole and energy using a gradient expansion of T^JEd as a 
framework. Hybrid models ,tlrll3 which have been partic- 
ularly successful in achieving quantitative improvements 
of binding and total energies in molecules, may also be 
motivated as approximations to the coupling constant in- 
tegration of the adiabatic connection formulaES 

In developing many of these approaches, many-body 
calculations of the exchange-correlation hole and related 
quantities in real systems, particularly in atoms and sim- 
ple molecules, have played a useful role, helping to con- 
firm assumptions of approximat£LfpK±ionals or point out 
specific areas for improvement .EZfEflO There have been 
few such studies for real solids, however 

A second fruitful decomposition of E xc has been to 
separate it into exchange and correlation components. 
The exchange energy, which contains the correlation due 
to Fermi statistics that occurs for the noninteracting or 
A = limit, may be explicitly written in terms of the oc- 
cupied single-particle orbitals obtained in this limit, and 
thus it and its functional derivative with density can in 
principle be obtained exactly and self consistently. Ap- 
plications along this line .-include the optimized effective 
potential (OEP) methodjij an implementation of the ex- 
act exchange potential which has been limited in prac- 
tice to free atoms and other simple geometries , and the 
method of Krieger, Li and Iafrate (KLI)E3 which pro- 
vides an approximate evaluation of the exchange func- 
tional by means of a simplified single-particle Green's 
function. Recently, an exact-exchange formalism appli-. 
cable to solids has been developed by Stadele et al£3'E£l 
These methods are of particular interest because they are 
self-interaction free and are promising casdidates for the 
quantitative treatment of excited states Jl3't3 

On the other hand, these advances expose an impor- 
tant problem: many traditional methods such as LDA 
rely for their success on the close relation between ex- 
change and correlation, in which the correlation hole 
tends to correct for, much of the nonlocal behavior in 
the exchange hole.E3 As a result, the error for the sum 
of exchange and correlation energies is typically signif- 



icantly smaller than for either separately! It is of in- 
terest, therefore, to study how the correlation hole alone 
behaves in real materials and how current functionals for 
correlation fare with respect to accurate many-body cal- 
culations. Or from the complementary point of view, it 
becomes necessary to consider to what extent modeling 
orbital effects in the correlation hola. become important 
once they are included in exchange O 

Variational Quantum Monte Carlo (VMC) calculations 
of the coupling-constant pair correlatiooJuaction for the 
Si crystal have been reported previously,Ej'E!J including an 
analysis in depth of the exchange-correlation energy den- 
sity obtained by the adiabatic connection theory. In this 
work we discuss in more detail the exchange-correlation 
hole, its decomposition into exchange and correlation and 
the relation between the angle-averaged hole and the en- 
ergy density. We present analysis along two lines: first 
a detailed comparison of the VMC exchange-correlation 
hole and that of several DFT models as a function of 
position in the crystal, in order to search for systematic 
trends in how these functionals describe the exchange- 
correlation hole of a prototypical semiconductor material. 
Secondly, we look at how the exchange and correlation 
holes separately behave in the bonding region of Si, and 
observe how orbital effects in the exchange generate re- 
lated effects in the correlation hole. 

The paper is organized as follows: Sec. [n] provides 
theoretical background on the exchange-correlation hole, 
and the various models used to describe it. We briefly 
discuss the method used to obtain computational results 
in Sec. 



Ill 



Results for spherically averaged holes are 
presented in Sec. [V, and a detailed analysis of the cor- 
relation and exchange holes in Sec. |y|. In Sec. VI we 
discuss issues raised by our data analysis concerning the 
improved treatment of nonlocal density information in 
DFT's, as well as orbital effects in the Si correlation and 
exchange holes that may be of relevance in improving 
upon exact ex chan ge methods. Our conclusions are pre- 
sented in Sec. VII. 



II. THE EXCHANGE-CORRELATION HOLE IN 
DENSITY FUNCTIONAL THEORY 

A. Basic theory 

We consider a family of systems parameterized by a 
coupling-constant A 



n-rj| 



(1) 



where V ex t is the external potential and V\ a potential 
added to keep the ground-state density invariant. The 
units here and elsewhere in the paper are in Hartree 



2 



REVTgX 3.0 



Cancio et al. Preprint, 2001 



atomic units unless otherwise indicated. Two limits of in- 
terest include the noninteracting system, A = 0, in which 
H reduces to the standard Kohn-Sham equation, with 
the adjustable potential equal to the Kohn-Sham poten- 
tial. The second is at A= 1, where V\ vanishes and one 
recovers the original fully interacting many-body system. 
The ground-state wave function is a Slater determinant 
of single-particle orbitals in the first case and the true 
many-body ground state in the second. 

The exchange-correlation hole is defined for a given 
value of A as the change in the ground-state expectation 
of the density at one point in the system r + u given the 
observation of an electron at some other point r: 

n xc (r,r + u;A) = 

7~t ( Y 6 ( T ~ S ( T + u ~ r j) ) a ~ n ( r + u )- 

Ti[Y ) 

(2) 

Here Q\ indicates an average over the ground-state wave 
function of H(X) and n(r) — {J2i^( r ~ r i))A is the A- 
invariant ground-state density. A related quantity, the 
pair correlation function g(r,r + u), is a measure of 
the pair density relative to that expected for uncorre- 
cted electrons with the same density distribution. The 
exchange-correlation hole is expressed in terms of g as 

n xc (v, r + u; A) = n(r + u) [ g(r, r + u, A) - 1 ]. (3) 

The importance of the exchange-correlatipti|Jjiole to 
density functional theory lies in its connectior£3'E3 to the 
exchange-correlation energy E xc : 



d 3 r n(r) 



dV 



dX 



c(r,r',A) 



(4) 



E xc includes all the contributions to the energy due 
to correlations, that is, beyond the noninteracting and 
Hartree energy. The exchange-correlation energy E xc [n] 
is analyzed in density functional theory in terms of 
the exchange-correlation energy density e xc (r) and the 
exchange-correlation energy per particle e xc (r) 



E xc = / d 3 r n{r)e xc (r) = / d 3 r e xc (r). 



(5) 



Relating this expression to that for E xc in terms of the 
exchange-correlation hole n xc one has 



e zc(r) = i J d 3 u 



1 f, 3 n xc (r,r + u) 



(6) 



where n xc (r, r + u) is the coupling-constant integrated 
hole, J Q dX n xc (r, r + u, A). The exchange-correlation en- 
ergy per particle thus has the natural interpretation as 
the interaction energy of the particle with its exchange- 
correlation hole. 

The exchange-correlation hole and energy are fre- 
quently decomposed according to separate exchange and 



correlation components. The exchange hole n x is that of 
the noninteracting or A = ground state [Eq. (||)] and 
contains correlations due solely to Fermi statistics. The 
remainder, n c , of the contributions to the total Tt X Q are 
those induced by introducing pair correlations into the 
ground-state wave function, representing density fluctu- 
ations caused by the Coulomb interaction. The two com- 
ponents satisfy important sum rules which reflect the 
conservation of particle number: 



J d 3 u n c (r,r + u) = 0. 



(7) 
(8) 



A key to the practical exploitation of the relation be- 
tween n xc and E xc is that, since the Coulomb interaction 
depends only on interparticle distance, much of the in- 
formation contained in n xc is not needed and may be 
ignored in the development of models. In particular, to 
obtain the energy per particle e xc (r) at a given point r 
the angular variation in n xc may be integrated out to 
leave a function depending only on the distance u from 



(n xc (r,u)) 



(9) 



with e xc (r) in Eq. |^ given in terms of the integral over u 
of the weighted angle-averaged hole, 2iru(n xc (r, u)). The 
DFT models we study in this paper can be viewed as 
models of this angle-averaged hole. 



B. Models for the exchange-correlation hole 

The LDA model for n xc may be obtained by replac- 
ing the true hole about r with that of the homogeneous 
electron gas at the local density, n(r): 



(n^ A (r,u)) =n(r){/ e >,r s (r)]-l}. 



(10) 



The factor r s = (3 /Ann) 1 / 3 is the average distance be- 
tween electrons in the homogeneous gas of density n. It 
provides a useful measure of the correlation length for the 
exchange-correlation hole, as the width of the exchange 
hole scales linearly with r s and that of the correlation 
hole shows similar behavior for the range of—densities 
1.4 < r s < 4 A valid for the Si crystal valence.o 

The WD Ad attempts to account for the inhomogeneity 
of the density in realistic systems by incorporating it in 
the density prefactor of n xc 

{nZ DA {*,u)) = J ^ n (r + u){g moda [u,f s (r)]-l}. 

(11) 

This form of the exchange-correlation hole treats the hole 
at r + u as the modification of the density at that point 
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rather than that at the center of the hole. This is a prop- 
erty of the exact hole that leads to a significant departure 
from the homogeneous case when the density varies con- 
siderably within the effective range of the hole. The ma- 
jor approximation made is the use of a scaling form for 
the pair correlation function where for each r an averaged 
length scale r s (r) is determined to select g from a fam- 
ily of suitably chosen pair correlation functions g(u,r s ). 
This length scale is fixed at each point fixed by the par- 
ticle sum rule 



(r + u){g m ° del {u,f s (r)) -1} = -1. (12) 



A natural choice for g(u, r s ) in an extended system is that 
of the homogeneous electron gas; this form will be con- 
sidered in this paper. The WDA scaling approximation 
has the ability to fit certain properties of the true n xc 
unattainable by the LDA; in particular it significantly 
improves on the description of the exchange-correlation 
energy per particle of atoms when a single electron is 
removed, obtaining the correct limiting value of 1/2R 
versus distance from the atom RE 

The ADAE3 is another method to incorporate into a 
model of the exchange-correlation hole information about 
the nonlocal density variation in the vicinity of the hole. 
It preserves the homogeneous gas model for n xc but de- 
termines the r s factor at which it is evaluated from an 
average of the density in the neighborhood of the refer- 
ence point: 



(n^(r,«)> =n(r){g h ^[u,f s (v)}-l} (13) 



with 



n(r) = / d 3 r'n(r') w(\r — r'|, r s ). 



(14) 



Here f s = (3/47rfi) 1 / 3 and w is a normalized weighting 
factor. The weight w is chosen to obtain the correct E xc 
in the linear response limit and has an effective range 
set by f s . The averaging procedure can in this sense 
be understood as the selecting of a physically reasonable 
correlation length for the hole from an average of the 
density over the physical range of the hole, particularly 
of value in cases where the density varies rapidly within 
this range. Unlike the WDA, the ADA does not try to 
provide a correction to the shape of the n xc . 

Note that with the use of a radially isotropic g, only a 
spherically averaged density enters into the definition of, 
and sum rule for, the WDA exchange-correlation hole: 



(nZ DA (r,u)) = (n(r,u)){g he °[u,r s (r)} - 1}, (15) 



and 



(n(r,u)} = { ^^-n(r + u) 

J 47T 



(16) 



This averaging reduces the information about the nonlo- 
cal density actually used by these methods from the full 
complexity of the density variation about each reference 
point to a simpler and often slowly varying radial func- 
tion (n(r,u)). Knowledge of this function at each point 
in space, along with the model used for g in the WDA and 
w in the ADA, specifies E xc within each approximation. 

Another line of attack in improving the exchange- 
correlation hole and energy, the GGA, has been to de- 
sign corrections to the LDA hole using information about 
system inhomogeneity contained in the local .density gra- 
dient. One implementation of this approacip has led to 
the construction of a GGA model for n^E^I In this and 
other GGA models the measure of inhomogeneity is the 
absolute value of the gradient, |Vn(r)|/n(r). Second- 
order effects proportional to the Laplacian of the density, 
V 2 n(r)/n(r), are mapped to a gradient form by an in- 
tegration by parts. This transformation has no effect on 
system-averaged quantities, namely E xc and the average 
over r of n xc . However, the local quantities n xc (r, r + u) 
and e xc (r) in the GGA are as a result no longer defined 
in terms of the adiabatic formulae of Eqs. Q aad || and 
not directly comparable with our VMC data.Ea 



III. SYSTEM AND COMPUTATIONAL DETAILS 

A description of the computational method used in 
obtaining numerical many-body expectations of the pair 
correlation function and related-quantities in Si has been 
presented in a previous paper .Ell We present here a short 
summary of the details relevant to our discussion. 

For V ex t in Eq. (|l|) we use a norm conserving nonlocal 
LDA pscudopotcntial which replaces the electrons in the 
atom core. The system and expectations involve valence 
electrons and their correlations. A simulation cell con- 
sisting of a 3 x 3 x 3 lattice of fee primitive cells of the di- 
amond lattice containing 216 valence electrons is used in 
obtaining expectations. Finite size effects in the electron- 
electron interaction were taken into account by replacing 
the Coulomb interaction in Eq. (Q) with a truncated po- 
tential into which the long-range Hartree interaction is 
incorporated. This form has proved to be successful in 
reducing finite size errors in the Coulomb enexgy as com- 
pared to the conventional Ewald interaction.Ej 

To describe the ground state for a given value of cou- 
pling constant, a Slater- Jastrow type wave function was 
used, consisting of a product of Slater determinants mul- 
tiplied by a trial many-body correlation factor: 

V^ A = cxp (-J2u x (\r i -r j \)+J2x X (^)) D]D[. 



(17) 



A similar result holds for the ADA, as the density averag- 
ing function w is a function only of interparticle distance. 



The orbitals used in the Slater determinants for up and 
down spins D a are determined from an LDA calcula- 
tion and include all the valence-band Bloch orbitals of 
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the primitive cell periodic on the supercell. The Jastrow 
factor introduces variationally explicit electron-electron 
correlations with one-body corrections to optimize the 
density. It is described in detail in the paper by Hood et 
alfl A variational ground-state wave function optimized 
over 22 variational parameters was obtained for five dif- 
ferent values of the coupling-constant, and recovered 88% 
of the correlation energy as compared to nearly exact dif- 
fusion Monte Carlo calculations for the same system. 

The pair correlation function g(r, r') [Eq.(||)] was ob- 
tained from the ground-state wave function for each A in 
order to obtain coupling-constant integrated exchange- 
correlation holes and energies. The pair correlation func- 
tion was expanded in symmetrized plane-waves, taking 
full advantage of the symmetries of the pair correlation 
function in the crystal to reduce the basis set to a reason- 
able number. Expectations for the coefficients for sym- 
metrized plane-waves were evaluated statistically using 
the Monte Carlo method. The resulting data are limited 
by plane-wave cutoff and statistical errors which can be 
checked by comparing the numerically calculated and ex- 
act values of the exchange hole. The resulting variation 
in g x was between 1% and 6%. The error due to basis 
set cutoff is most noticeable near the atom core where, 
despite the pseudopotential treatment, the rapid varia- 
tion of the density and the valence orbitals leads to large 
wave-numbers in the exchange hole. 



IV. ANGLE-AVERAGED 
EXCHANGE-CORRELATION HOLE 

Figure [l] shows as a function of interelectron distance u 
the angle-averaged correlation, exchange and exchange- 
correlation holes, weighted by 2nu, for the VMC calcula- 
tion of Hood et aLLEEi and several DFT models at sev- 
eral representative points in the Si crystal. The pair cor- 
relation function of the homogeneous electron gas used 
to generate the DFT holes is obtained from an accurate 
model form.EZI The distance weighting factor 2ttu is cho- 
sen so that the integral of the plotted function yields the 
exchange-correlation energy per particle at the reference 
point. These values are listed separately in Table |[ 

In order to give an idea of the environment surrounding 
each reference point and to compare with the WDA and 
ADA forms for n xc , the angle-averaged density (n(r, u)), 
the average density at a distance u from the reference 
point r, is shown at the top of each plot. The density 
is to be compared to the average over the unit cell of 
0.0296 a.u. The angle-averaged density for each of the 
four reference points shown indicate four different envi- 
ronments typical of the pseudopotential model Si crystal. 
In the first, (a), the reference electron is at a density con- 
siderably higher than the unit-cell average, and the den- 
sity decreases significantly over the length-scale of the 
VMC exchange-correlation hole. The antibond position 
(b) is typical of intermediate density points where the lo- 



cal density is close to the unit-cell average and the angle- 
averaged density in the vicinity of the reference point 
shows only modest variations over long distances. The 
atom center (c) is typical of the "pathological" situa- 
tion inside the pseudopotential core. Here the density 
increases by an order of magnitude, from roughly 10% 
to 200% of the unit-cell average, over roughly one-third 
of the Si bond length of 4.44 a.u., before slowly decay- 
ing back to the average. Such an environment of dra- 
matic density change over the length scale of the hole 
indicates the likely failure of the LDA in this region, 
which has in fact been observed— in the comparison of 
LDA and VMC energy densities .Ej Finally, the intersti- 
tial point (d) shows the generic features at low density 
outside both valence and core. Here the angle-averaged 
density also changes considerably from the local density 
but at a much more gradual rate. 

In the high-density case of the bond center, Fig. [l](a), 
the exchange-correlation hole closely follows that deter- 
mined by the local density, with r s = 1.397 a.u. The aver- 
age density within the effective range of the hole is lower 
than the local density, leading to a slightly larger f s for 
the ADA and an overall better fit of the hole. The WDA 
hole [Eq. (|l5|)] shows the effect of two competing factors. 
The value of f s used in the scaled pair correlation func- 
tion g heg is significantly larger than the local r s . The 
resulting pair correlation function is wider and deeper 
than that of the LDA. However, the density prefactor 
(n(u)) in Eq. ( p"5| ) drops rather quickly and "quenches" 
the exchange-correlation hole at larger u. Eventually, as 
the density settles to its asymptotic value of the unit cell 
average, the hole has a long-range tail again dominated 
by the shape of the pair correlation function. The net 
result is a considerable overcorrection of the LDA model 
with respect to the VMC data. 

At the antibond point, the hole is wider and less deep 
than at the bond center, consistent with the variation 
of the homogeneous electron gas hole with respect to 
density. Neither form of nonlocal averaging results in 
a significant change from the LDA form of n xc since the 
angle-averaged density is almost constant in the vicin- 
ity of the reference electron. The exchange-correlation 
energy per particle, determined from the integral of the 
curves in Fig. |l], has roughly the same modest deviation 
from the VMC data in each shown in Table Q. 

This leads to a significant difference in the energy density 
because the density is still large at this point. 

The holes at the bond center and antibond point 
are fairly representative of those at other points of 
intermediate-to-high density in the Si crystal. The cor- 
rections introduced by the nonlocal DFT's are quite small 
because of the relatively small variation in the angle- 
averaged density within the effective range of the hole. 
These high density corrections to n xc and e xc have a sen- 
sitive dependence on position in the crystal, which ap- 
pears to lack an obvious relation to the positional trends 
in the VMC data. As a result, plots of the deviation of 
the ADA and WDA energy density from the VMC data 
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FIG. 1. Distance- weighted angle-averaged exchange-correlation holes at various points in the Si crystal: (a) bond center, 
(b) antibond position, (c) atom center and (d) tetrahedral interstitial site in the (110) plane. Correlation, exchange and 
exchange-correlation components for VMC data (thick solid lines) and several models are shown, with the exchange set off by 
-0.05 a.u. and exchange-correlation by -0.15 a.u. on the vertical axis for clarity. Solid line on the top half of each subplot shows 
the angle-averaged density langlen(r, u)) as function of distance from reference electron. All quantities plotted are in atomic 
units. 
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TABLE I. Exchange-correlation energy per particle e xc (r) in hartees at various points in the Si crystal discussed in the 
text. Label VMC refers to raw data, VMC-I includes a correction to obtain the correct particle sum rule, VMC-II includes a 
correction for the plane-wave cutoff (see text) and corresponds to the integral over the angle-averaged holes shown in Fig. hi. 
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feature quite complex "dapple" patterns that remain af- 
ter jtlj^-jworst errors of the LDA model are smoothed 

For the two low-density cases, however, the nonlocal 
density averaging techniques have a considerable impact 
on the shape and quality of the exchange-correlation hole. 
In these cases the shape of hole is dominated by the rapid 
increase of the density with u. 

At the atom center, Fig|l|(c), the rapid change in den- 
sity leads to a far more compact exchange-correlation 
hole than that predicted by the LDA at the low local den- 
sity (r s = 4.12 a.u.). The minimum in the weighted hole, 
where it contributes the most to the exchange-correlation 
energy, is 1.6 a.u. from the reference point, less than half 
that of the LDA hole. Roughly 90% of the total contri- 
bution to e xc (r) comes from within a radius of 3.0 a.u., 
significantly less than a bond length (4.44 a.u.). In this 
case, both ADA and WDA do significantly better at de- 
termining the overall depth and length scale of the hole 
over the energetically important region of 0.5 to 3.5 a.u. 
In addition, WDA clearly matches the shape of the VMC 
hole quite closely. As a result, the ADA and WDA obtain 
fairly accurate estimates of the energy of the exchange- 
correlation hole as shown in Table Q while LDA obtains 
only half the VMC value. 

The interstitial point, Fig. |l|(d), has low-density fea- 
tures qualitatively similar to the atom center but, con- 
sistent with a more gradual change in density, is consid- 
erably shallower and more spread out. The WDA has a 
particularly good overall agreement with the VMC hole 
on a point-by-point basis (although the ADA has a total 
energy slightly closer to the VMC value). However, possi- 
bly due to the more gradual change in density within the 
vicinity of the hole, the LDA does a much better job of 
matching the energy of the hole than near the atom cen- 
ter. The large differences in the shape of the hole largely 
cancel out in the integral so that the energy of the hole 
is only 20% smaller than the VMC value. The error in 
the LDA energy density, which is obtained by weighting 



the energy of the LDA exchange-correlation hole by the 
density, is thus fairly insignificant at this point :cj in con- 
trast the LDA energy density in the atom core suffers 
from considerably larger deviations with the VMC. 

Figure [l| also shows results for angle-averaged correla- 
tion and exchange holes at each point. In general, the 
exchange hole is the dominant contribution to the total 
hole. The correlation hole provides a correction to the ex- 
change hole, reducing electron density near the reference 
electron and enhancing it towards the outside of the ex- 
change hole, thereby making the total hole slightly deeper 
and more compact. At high densities this correlation hole 
is weaker and shorter-ranged than the LDA hole, a re- 
sult in agreement with the prediction of gradient based 
models£3 This trend however does not carry over to low 
densities, particularly in the atomic core [Fig. |l](c)]. The 
correlation hole here, while clearly shorter-ranged than 
the LDA hole, is also significantly larger and makes a 
much larger relative contribution to e xc (r). The WDA 
and ADA both do moderately well in predicting the mag- 
nitude and length scale of the correlation hole in this case. 

In general the nonlocal density averaging methods fair 
less well at predicting correlation or exchange alone than 
the combination of the two. For the WDA in particular, 
the sum rules of the correlation hole and the exchange 
hole alone are not satisfied by satisfying the sum rule for 
exchange-correlation [Eq. (12)], with the correlation hole 
influenced more by the short-range behavior of the angle- 
averaged density and exchange the long-range behavior. 
This may contribute to the larger errors seen in n c and 
n x than for n xc observable in the cases with significant 
density variation. 

Some sense of the general trends resulting from the 
different strategies for choosing the pair correlation func- 
tion can be obtained from plotting the solution for the 
average interelectron distance f s used by the nonlocal 
methods versus the r s obtained from the local density. 
The result for sample points at various densities through- 
out the unit cell is shown in Fig. ||. In the case of the 
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FIG. 2. Correlation lengths of WDA and ADA ex- 
change-correlation holes. The correlation length for the equiv- 
alent homogeneous electron gas hole, f a = (3/47rn(r))~ 1/ ' 3 , 
obtained from the nonlocal DFT models discussed in the pa- 
per is plotted versus the correlation length r s derived from the 
local density for holes at representative points throughout the 
unit cell. Dotted line gives the LDA approximation, dashed 
line the r s factor obtained from using the unit cell average. 



ADA, for r s lower than the unit-cell average of 2.01 a.u., 
f s hews closely to the local density value. It deviates 
from it at lower densities, either dramatically at the atom 
center (the one point below the average r s curve), or 
more gradually as one heads towards the interstitial point 
(the other low-density points plotted) . The WDA values 
in contrast are all grouped about the unit-cell average 
value, regardless of location and density. In other words, 
the variation in the WDA exchange-correlation hole is 
primarily derived from variation in the nonlocal angle- 
averaged density, and very little from the variation of g. 
These results reflect the character of the criteria used to 
determine the average density in each case. The weight- 
ing function of the ADA is an oscillatory function in real 
space, with a peak contribution at roughly 0.5 r s from the 
reference electron, so that a significant effect is observed 
only where the density variation about the electron is 
sufficiently rapid. In contrast, the weight implied by the 
WDA sum rule condition, u 2 [l — g heg (u,r s )], is peaked 
near r s and has significant contributions out to 2 r s , a 
distribution broad enough to to pick up the average den- 
sity at almost every point in the crystal. 



V. EXCHANGE AND CORRELATION HOLES 

As discussed earlier, the trends in the comparison of 
the DFT and VMC models for e xc (r) are complex and 



difficult to characterize. Some insight into what is go- 
ing on can be obtained by decomposing the exchange- 
correlation hole into exchange and correlation compo- 
nents and particularly by considering the response not 
only as a function of distance from the reference electron 
but including the complete information, n xc (r, r + u), of 
the second electron's position in the crystal. 

The exchange hole is of use particularly in that it can 
be derived exactly and thus provide an unambiguous test 
of density functionals. Further, it is often the dominant 
part of the total exchange-correlation hole. It may be 
formally obtained as the exchange-correlation hole as- 
sociated with the Slater determinant ground-state wave 
function of the noninteracting (A = 0) limit of Eq. [j]. By 
applying this wave function in Eq. one has 



n x (r,r + u) 



n(r) 



E 



^Va(r)C(r + u) 



(18) 



where o is a spin index and N a the number of electrons 
with spin a. In practice n x is obtained from LDA va- 
lence orbitals periodic on the 3x3x3 simulation cell 
used in our VMC calculation. These produce a density 
that is indistinguishable from the VMC density within 
the statistical limitations of the Monte Carlo sampling. 
At u = the exchange hole reduces to — n(r)/2, which 
reflects the Pauli exclusion prohibiting the occupation of 
the same point in space by two particles of the same spin. 
In addition, the exchange hole is strictly negative, as may 
easily be seen in the case of Si where the orbitals may be 
taken as real. 

In Fig. H we show results for the exchange hole and 
angle-averaged hole weighted by distance for three po- 
sitions near the bond center of the Si crystal. In the 
bottom half of the figure is a contour plot of n x (r, r + u): 
the exchange hole at r + u given a reference electron at 
the position r. The plots cut through a chain of Sili- 
con atoms in the (110) plane of the crystal; the reference 
electron's position is shown as a white dot. The shad- 
ing represents a change in value ranging from just less 
than zero for the white region to just above the absolute 
minimum of minus one-half the peak density of 0.087 
a.u. in the black regions. The thick solid curve in the 
upper plots show the result of spherically averaging the 
exchange hole to obtain a function of distance from the 
reference point. The same weighting factor of 2ttu as that 
of Fig. |l| is used. This result is compared to the various 
DFT models discussed previously. 

In the case of the reference electron at bond center, the 
exchange hole is focused in the bonding region immedi- 
ately surrounding the reference electron, with the bulk 
of the hole contained in the region between the atoms 
on either side of the bond, and a negligible weight on 
the nearest neighbor bonds. The electron density has 
been reduced by roughly 45% within the bond, i.e., in a 
diamond-shaped region extending 1.7 a.u. along the bond 
axis in either direction and 1.2 a.u. perpendicular to it. 
This accounts for 90% of the density of electrons with 
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FIG. 3. Exchange hole near bond center: (a) off bond, perpendicular to bond axis, (b) bond center, (c) off bond center along 
[111] axis. Top half shows angle- averaged hole {n x (u, r)) weighted by 2ttu versus distance from reference electron u for VMC 
data and several models. Bottom half shows contour plot of n x (r, r + u) in (110) plane, with location of reference electron as 
white dot. Gray scale is in increments of 0.005 a.u. with white region between -0.005 and 0.000 a.u. 



the same spin as the reference one. As one moves the 
reference electron either along the bond axis (c) or per- 
pendicular (a) to it, the hole tends to remain focused on 
the nearby bond center, despite a small shift in weight 
towards the reference point. As before, more than 90% of 
the same-spin density is removed in the bond. This stiff- 
ness or insensitivity to the position of the reference elec- 
tron is gradually lost as the reference electron is moved, 
say, to the atom center or towards the antibond point 
equidistant from three bonds. In the latter case, the ex- 
change hole is again centered on the reference point and 
resembles a sp 3 -hybrid atomic orbital. 

The effect of the stiffness of the exchange hole on its an- 
gle average can be seen in the difference between the LDA 
and exact cases (thin and thick solid lines). In the bond 
center, Fig. |^(b), the LDA result closely matches the ex- 
act value. At position (a) off the bond axis, the bond 
center is oriented tangentially with respect to an angle 
average about the reference point and the exchange hole 
has become deeper and narrower than that obtained by 
the LDA. In contrast, in case (c), where the bond cen- 
ter is oriented radially out from the reference electron, 



and the electron is near to an atomic core, the hole has 
become shallower and wider than in the LDA. Although 
the two positions have similar densities and density gra- 
dients, the resulting deviations from the LDA are qualita- 
tively different and point to a genuinely nonlocal behavior 
in the exchange hole. Nonetheless, neither the WDA nor 
the ADA are particularly sensitive to the position of the 
bond center with respect to the reference electron and 
provide no systematic improvement (or disimprovemcnt) 
over the LDA fit. 

The correlation holes corresponding to each case in 
Fig. H are shown in Fig. |]. Note the order of magnitude 
smaller range of density changes involved. The corre- 
lation hole n c takes on both negative values in regions 
where electrons are repulsed by the Coulomb potential 
of the reference electron and positive values representing 
regions where electron density has increased. The net 
electron density change is zero. 

At the bond center, Fig. [|(b), the correlation hole con- 
tour plot can be roughly separated into two regimes: a 
deep and narrow well in the vicinity of the reference elec- 
tron and a shallow longer-ranged response that may be 
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FIG. 4. Correlation hole near bond center: (a) off bond, perpendicular to bond axis, (b) bond center, (c) off bond center 
along [111] axis. Top half shows angle-averaged hole weighted by 2nu versus distance from reference electron u for VMC data 
and several models . Bottom half shows contour plot of n c in (110) plane, with location of reference electron as white dot. 
Gray scale is in increments of 0.001 a.u. with white region between 0.001 and 0.002 a.u. 



described as a polarization of the nearest neighbor bonds. 
In the region of the pseudopotential core the correlation 
response is suppressed. The resulting angle-averaged hole 
has a complex and unusual structure as a function of dis- 
tance. In comparison, the LDA and other DFT holes 
show a relative lack of structure and greatly overestimate 
the magnitude of the correlation energy per particle at 
bond center. 

As the position of the reference electron is moved either 
off the bond axis (a) or along the axis (c) , the correlation 
hole demonstrates a marked sensitivity to its position 
relative to the bond center. In either case the minimum 
is slightly offset in the direction of steepest increase in 
density. This is consistent with g c being nearly isotropic 
near the electron as expected from the cusp condition. 
Outside this short-range region the hole undergoes strong 
distortions in shape from that of the hole about the bond 
center. These longer wavelength features, as in the bond- 
center case, are well characterized in terms of the polar- 
ization of the nearby bonds. Bond polarization is particu- 
larly noticeable for the electron fixed on the bond axis (c) 
where there is a large shift of the electron density from 



the atom nearest the electron to the other side of the 
bond, with a peak on the bond axis opposite the refer- 
ence electron. Nearest neighbor bonds are polarized in a 
similar fashion. In (a), the nearby bonds are polarized 
with density shifted from the nearer side of the bond with 
respect to the electron position to the farther, with the 
node of the polarization approximately along the bond 
axis. Consideration of the corresponding pair correlation 
function shows that these features in the correlation hole 
constitute a significant departure from an isotropic form 
of pair correlation function such as is assumed by the 
DFT models discussed in this paper. 

Despite the reduction of information engendered by 
averaging over angle, the dramatic changes in the shape 
of n c with position show up in the angle-averaged hole, 
giving rise to nonlocal corrections to the LDA hole. The 
LDA model depicts a trend in n c to a broader and shal- 
lower hole as the density decreases from case (b) to (c) 
to (a) in the upper half of Fig. g. The VMC hole in 
case (a) retains the character of the bond center hole, 
being narrower and weaker than that predicted by the 
LDA. The overall disagreement with LDA is less pro- 
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nounced and some of the detailed structure is lost. As 
one moves along the bond axis to case (c), however, a sig- 
nificantly different trend occurs: the hole in the region of 
density reduction becomes deeper and broader relative 
to the LDA, and in particular now matches the width 
of the LDA hole. In addition to the stronger response 
in this region, a large peak appears at the hole edge in 
keeping with the zero sum rule of the correlation hole. 
This peak, at about 3.5 a.u., correlates with the position 
of electron density peaks on the nearest neighbor bonds 
in the contour plot. As with the case of exchange, these 
trends in n c poorly correlate to the changes in density 
or density gradient, and show a sensitivity to the larger 
nonlocal crystal environment. 

Though the exchange and correlation holes both show 
significant and nonlocal deviations from the isotropic 
form typical of the homogeneous electron gas, these de- 
viations are highly correlated with each other and thus 
tend to cancel in their sum. As the exchange hole tends 
to stay on the bond center, the bond polarization in the 
correlation hole is oriented so as to shift the center of the 
exchange-correlation hole towards the reference electron. 
A simple picture of this feature in the correlation hole is 
that it represents a response to the anisotropic distribu- 
tion of the exchange hole by which some of the deviations 
of the exchange hole from an isotropic form are screened 
out in the correlation hole. This behavior may help to 
explain the relative success of the LDA in describing the 
exchange-correlation hole. 

Neither the LDA nor the nonlocal DFT's model with 
much fidelity the variation in the VMC correlation hole, 
aside from the general broadening of the hole with lower 
density. In particular the corrections of the WDA lead 
to a worse fit of both the correlation hole and its en- 
ergy. The differences in exchange and correlation be- 
tween WDA and VMC do tend to balance each other 
but not consistently: for example they cancel nicely in (a) 
and (b) but add in case (c). The net effect of these trends 
on the exchange-correlation energy is shown in Table ||, 
where the total E xc , E x and E c for several methods are 
compared. In addition to the methods discussed in detail 
in this paper, we present results from the PW91 version 
of the GGA land those of a diffusion Monte Carlo (DMC) 
calculation.ESa The DMC method removes nearly all of 
the variational bias in the VMC correlation energy at a 
considerably larger computational cost. With the excep- 
tion of the GGA, none of the density functional meth- 
ods obtain a good estimate for E x or E c as compared to 
DMC, with errors of roughly 1 to 2 eV per atom. Esti- 
mates for E xc are much closer, especially for the WDA. 



VI. DISCUSSION 
A. Limits of density averaging 

The clear success of density averaging occurs at 
low densities where ADA and WDA obtain excellent 
exchange-correlation energy densities. WDA in par- 
ticular predicts the shape of the exchange-correlation 
hole with exceptional fidelity, even in the extreme sit- 
uation inside the atomic core. At high density, the small 
level of variation in the nonlocal angle-averaged density 
[Eq. (|l^)] limits the effects of density averaging to subtle 
alterations of the hole which do not provide a systematic 
improvement with respect to the LDA. More importantly, 
the discrepancy between the density-averaged holes and 
the VMC hole (much less the exact hole) are difficult to 
correlate with any known quantity. 

However, density averaging must lead to a global and 
systematic improvement over the unit cell or else the 
quality of the total result may be disappointing. This 
is particularly the case of the ADA, which returns val- 
ues very close to the VMC values at low density and 
slightly underestimates the magnitude of E xc at high 
density. Unfortunately, once the variational bias of the 
VMC energy is removed by a DMC calculation as shown 
in Table || this result proves to be a significant disim- 
provement with respect to the LDA. Moreover, although 
the different prescriptions for the WDA and LDA holes 
lead to extreme differences in the degree of nonlocality 
they incorporate, both end up with quite similar, rea- 
sonable predictions for the total E xc . Apparently, the 
large error at low densities of replacing n(r + u) with the 
local density n(r) in the definition of the LDA exchange- 
correlation hole is offset by correspondingly large changes 
in g, which are suppressed by density averaging in the 
WDA. 

It is interesting to compare the density averaging ap- 
proach to the gradient expansion used as the basis for the 
GGA. A notable result of the PW91 version of the GGA 
is that it not only provides a significant improvement 
of the LDA exchange-correlation energy, but of the ex- 
change and correlation energies as well (Table ||). Thus, 
it alone among the functionals we have studied shows 
promise to be a reasonable candidate as a correlation- 
only density functional for this material. In comparison, 
the WDA, although it returns a very good E xc , has the 
worst estimate for E c . However, the local density and its 
Laplacian, the formal equivalent of the information about 
the surroundings of an electron used in obtaining a GGA 
hole, constitutes but a small part of what is contained 
in the angle-averaged density {n(r, u)}. The limitation of 
the WDA is the restriction to a pair correlation function 
form with only one variable parameter, one to which the 
hole in our case proves largely insensitive. It would be in- 
teresting if more of the information contained in (n(r, u)) 
could be used as input into a more flexible model for g, 
particularly for correlation or exchange separately. Possi- 
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TABLE II. Exchange and correlation energies E x and E c in eV per atom for VMC and various density functional methods 
described in the text. The results of a diffusion Monte Carlo DMC calculation to remove the variational bias due to the use of 
a trial wave function is also shown. 





LDA 


ADA 


WDA 


GGA 


VMC 


DMC 


E x 


-27.66 


-27.56 


-27.37 


-29.10 


-29.15 


-29.15 


E c 


-5.09 


-5.11 


-5.63 


-3.93 


-3.58 


-4.08 


E xc 


-32.75 


-32.67 


-33.00 


-33.03 


-32.73±0.01 


-33.23±0.08 



bly, the rate at which the density changes over the length 
scale of the hole, a factor which, as illustrated in Fig. [j], 
is clearly important in determining the overall shape of 
the WDA hole, could be used to influence the form of g 
in analogy to the GGA. 

Furthermore, (n(r, u)) has truly nonlocal information 
not accessible to GGA's. A salient feature of the system 
in our study is the long-range limit of the angle-averaged 
density, which tends to a finite constant for a crystal 
and to zero for an atom or molecule. It is thus possi- 
ble that low density points in these two systems have 
similar local or semilocal environments (and thus the 
same GGA holes) but significantly different exchange- 
correlation holes because of different long-range bound- 
ary conditions. Because of the nonlocal character of the 
WDA, these boundary conditions have a strong influence 
on the WDA hole, causing the length scale f s to be nearly 
constant for the crystal but for the atom, to range from a 
finite value near the valence density maximum to infinity 
as the reference electron is removed from the atom. The 
WDA produces a reasonable fit to n xc at low density for 
both the Si atom and Si crystal, despite the markedly 
different f s factors necessary to satisfy the particle sum 
ruleo It would seem, then, that the WDA should be 
capable to match or surpass the quality of the GGA - 
perhaps given a more accurate and flexible model pair 
correlation function form. 



gradient corrections to the LDA on a system-averaged 
basis. It does in fact capture general trends in exchange 
and correlation as discussed briefly in the next section, 
but clearly not on a local or point-by-point basis. 

This approach is frequently justified by the observa- 
tion that the energy density is not in itself a necessary 
component of an accurate DFT, rather only the total 
E xc and its functional derivative with respect to the den- 
sity, the exchange-correlation potential.El Nonetheless, it 
is worth mentioning that the deviation of the LDA energy 
density from the VMC energy density closely follows the 
Laplacian of the density. For example, the minimum of 
the Laplacian, an hourglass-shaped region near the bond 
center, closely correlates with the region where we have 
found the largest positive deviation of the LDA energy 
density from the VMC valucE] Likewise a large nega- 
tive deviation in the LDA energy density occurs in the 
atom center where the Laplacian has a maximum, and a 
smaller negative deviation in the interstitial region where 
the Laplacian has a weak local maximum. This latter 
feature is clearly visible in the n xc and e^r) data of 



Sec. [V. We expect that a GGA model that fits both 



gradient and Laplacian terms could prove very useful in 
improving the energy density and thus also the exchange- 
correlation potential. 



C. Orbital effects in exchange and correlation 



B. The exchange-correlation hole and energy density 
in the GGA 

As discussed in Sec. ||, the GGA models gradient ex- 
pansion corrections to the LDA in terms of the gradient 
of the density alone, mapping corrections proportional to 
its Laplacian to a gradient squared term by an integration 
by parts in Eq. This transformation has significant im- 
plications for the GGA models of n xc . For example, each 
position in Fig. |l| is a critical point in the density, and 
as such characterizes the nature of n xc in its vicinity. At 
these points, the PBE model for n xc is indistinguishable 
from that of the LDA, a particularly serious error in the 
atom core region |l](c) , where effects of inhomogeneity are 
most apparent. At the same time, corrections to n xc from 
the atom core and other regions where the Laplacian is 
large are mapped to the n xc of regions where the gra- 
dient is large, producing misleading corrections in these 
regions as well. The GGA is rather designed to provide 



Insight into the comparison of our VMC results and 
those obtained from the various models derived from the 
homogeneous electron gas may be obtained from a con- 
sideration of the semiconductor environment. 

First of all, a distinguishing feature of systems such as 
molecules or semiconductors which have a finite energy 
gap is that the ground state may be described in terms of 
an exponentially decaying localized basis. In a periodic 
system these are the Wannier functions, which in Si can 
be defined by a unitary transformation from the four va- 
lence bands to four prhitals per unit cell localized about 
the bond centers r/:cilc3 

Wl{T -n)=J2 Unk,l4>nk{v)- (19) 
nk 

The symmetry of the crystal requires that each orbital 
W[ must be related to the others by a space group op- 
eration. Such orbitals are well localized on a given bond 
and describe well the bonding character of the valence 
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FIG. 5. VMC exchange and correlation energy per particle in the (110) plane of the Si crystal, relative to the LDA value. 
The contours show increments of 0.01 a.u., with the thicker contour showing where the VMC and LDA values are equal. White 
dots show the locations of the reference electrons about which exchange and correlation holes are shown in Figs. ^ and |^. 



electrons. The exchange hole is most clearly represented 
in terms of these local functions and deviates strongly 
from the homogeneous electron gas model when these 
functions have little overlap. The exchange hole in terms 
of a Wannier basis is 



n x (r,r + u) = 

--Y 

n(v) ^ 



^WMr- r/ )W7(r + u-r 7 ; 



(20) 



Using Wannier functions obtained from our_LDA orbitals 
and the method of Marzari and Vanderbilt pil we find that 
97% of the the on-top value n x (ri,ri) of the exchange 
hole at bond center is determined solely from the Wan- 
nier function localized on that bond site. Thus, the lack 
of sensitivity of the exchange hole with respect to refer- 
ence position seen in Fig. || near the bond center is likely 
a reflection of the domination of the exchange hole in 
this region by a single Wannier function. This situation 
is reminiscent of that of the H2 molecule in which the 
exchange hole is constructed from a single orbital and is 
totally insensitive to electron position. 

An exponentially decaying exchange hole typical of an 
insulator or finite system has the effect of lowering the 
exchange energy with respect to that of the LDA, owing 
to the more localized form of the hole as compared to 
that of the homogeneous electron gas; imposing a finite- 
ranged hole has been m ingredient in constructing suc- 
cessful GGA methodsEj In our case the actual E x is 1.5 
eV lower than the LDA result, in good agreement with 
the PW91 GGA result; likewise, the the exchange en- 
ergy per particle e x (r) is lower than the LDA prediction 
through much, but not all, of the unit cell as shown in 
Fig. | 



The semiconductor environment also affects the corre- 
lation hole in several ways. First, there is a finite en- 
ergy cost to correlate electrons as compared to the ho- 
mogeneous electron gas and therefore there should be a 
smaller correlation response and correlation energy than 
predicted in the LDA. This is seen in the VMC corre- 
lation energy density and energy per particle at almost 
all points in the unit cell, as shown in Fig. ||. Again 
this trend is consistent with the assumptions made in 
the PW91 version of the GGA but does not correlate 
with those made in the density averaging methods. 

Secondly, as the electronic structure of the covalent 
bond plays an important role in exchange (in the form of 
a bond-orbital- like exchange hole near the bond center), 
it should do so as well in correlation. The signature of 
the existence of the covalent bond in the correlation hole 
is bond polarization, a clear feature .of our VMC data. 
In the simplest bond-orbital picture^ this arises from 
the introduction into the noninteracting ground state of 
excited states in which a pair of electrons are excited 
from-bonding to antibonding orbitals on the same bond 
site.Ej These states contribute a "left-right" correlation 
similar to that of H2 and other diatomic molecules in 
which the probability density of the second electron is 
shifted along the bond axis to the opposite side of the 
bond as that occupied by the first electron. A similar 
excitation of nearest neighbor electrons causes van der 
Waals like correlation between bondso 

These effects are evident in the "bond-right" case (c) 
of Fig. H where the antibond orbital has a peak. At the 
same time, the correlation hole associated with this left- 
right polarization is negligible when the reference elec- 
tron is placed on the nodal plane of the antibond which 
passes through the bond center normal to the bond axis: 
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that is, the situation of cases (a) and (b) of Fig. [|. Con- 
sequently, the weakness of n c relative to the homoge- 
nous electron gas based DFT models in these two cases 
and the enhancement of the correlation hole on the bond 
axis have a plausible explanation in the differing contri- 
bution of the antibond correlation to the holes in each 
case.c3 Moreover, bond polarization contributes notice- 
ably to the functional variation of e c (r) in this bonding 
region. We find that the e c (r) obtained from VMC data 
has a maximum (that is, a minimum in the energy re- 
duction obtained by electron correlations) on a narrow 
ridge centered on the antibond node passing through the 
center of the bond with normal parallel to the bond axis. 
Although some of the variation in e c (r) in the bonding 
region can be accounted for within the LDA, the behavior 
is sufficiently dissimilar to the local density variation as 
to lead to the largest discrepancy in the unit cell between 
the VMC and LDA results, as shown in Fig. ||. 

The significance of the antibond correlation and bond 
polarization in general is that it provides a partial expla- 
nation of the cancellation of errors in the LDA exchange 
and correlation hole seen in their sum. The exchange- 
correlation hole of the homogeneous electron gas is es- 
sentially "dynamic" in nature, that is, responding to the 
position of the reference electron and largely insensitive 
to the details of electronic structure. The "nondynamic" 
or orbital-dependent features in the exchange hole, such 
as its insensitivity to particle position near bond center 
are a large potential source of error for the LDA and 
other models based on the homogeneous electron gas. 
The partial cancellation of these effects by a correspond- 
ing nondynamic feature of the correlation hole should 
then provide a combined hole much more amenable to 
the LDA and related density gradient corrections. 

We have made preliminary estimates of the antibond 
contribution to the correlation hole n c and the correla- 
tion energy per particle e c (r) in the bonding region of 
Si using perturbation theory and perturbative configura- 
tion interaction methodsEj and localized orbitals derived 
from pseudopotential plane- wave DFT orbitals. We re- 
produce the qualitative differences in e c (r) between the 
points studied in this paper, as well as a reasonable re- 
production of the bond polarization features of Fig. ||. 
Details of this calculation are to be reported in a later 
paper. 



D. Error analysis 



As discussed in Section III, the exchange-correlation 
hole suffers from fairly large plane-wave cutoff errors in 
the atomic core where the pseudopotential orbitals vary 
rapidly. In order to produce a reasonable correlation hole, 
particularly in the core, it was necessary to correct for 
this error by taking the n x and n xc expanded to the same 
cutoff. Then a best estimate for n xc was taken by adding 
the finite cutoff estimate of n c to the exact n Xl that is, 



an estimate of n xc obtained from correcting the exchange 
component of the Monte Carlo estimate: 

n xc = n a x p c prox +n x - n a x pprox . (21) 

The change in n xc is noticeable particularly for the atom 
center case Fig. |l](c) , where the plane- wave correction re- 
sults in a deeper and more localized n xc . The resulting 
e xc (r), shown in Table |, varies from the raw VMC data 
by roughly 8%, increasing the agreement with WD A and 
ADA at the expense of the LDA value. However, given 
the poor agreement of the LDA and VMC values, the un- 
certainty in the VMC value docs not significantly alter 
the assessment of a dramatic improvement by ADA and 
WDA in this region. A rather smaller change was no- 
ticeable for the angle-averaged hole for the "bond-right" 
position discussed in Figs ||] and [l| In this case, the dif- 
ference is noticeable on the scale of the energy differ- 
ences between the various theories, with the plane- wave 
correction favoring the ADA case. The correction was 
essentially negligible for the other points studied. 

A second source of error in the VMC data is due to the 
statistical nature of the Monte Carlo estimates of expec- 
tations. Statistical sampling leads to roughly equal errors 
in the expectations of each measured plane-wave com- 
ponent of the pair correlation function, and a homoge- 
nous background noise (up to the resolution of the plane- 
wave expansion) in the real space behavior of the func- 
tion. With the assumption of homogeneous and uncor- 
related background noise, the weighted angle-averaged 
holes shown should have an error that is roughly con- 
stant at large interparticle distances. We find this to be 
a reasonable description of the long-range fluctuations in 
this quantity, resulting in a very small error for all the 
cases we have studied. On the other hand, the energy and 
particle sum rule of the exchange-correlation hole, both 
integrals over a large volume, suffer more serious cumu- 
lative effects from background noise, particularly at large 
interparticle distances where n xc is essentially zero. As 
a result, the VMC particle sum rule as estimated from 
the plots in Fig. |] typically varies by 1 to 3% from the 
correct value. The values for the exchange-correlation 
energy per partkle_and energy density reported in our 
previous studiesEH have been made by adjusting the 
exchange-correlation hole at each A to obtain the correct 
particle sum rule, and numerical values at the various 
points considered in detail here are shown in Table Q. 
This correction turns out to be noticeable at high den- 
sity relative to the small differences between VMC and 
the various DFT models. 

A final source of error is variational bias in the trial 
ground-state wavefunction resulting from its limited vari- 
ational flexibility. The resulting discrepancy from the 
true ground state can be significant for correlation: the 
optimized VMC correlation energy is 14% higher than 
that estimated by the nearly exact DMC calculation. Ex- 
change, on the other hand, is relatively unaffected by 
variational bias, as the VMC density and the associated 
single-particle orbitals are expected to be very accurate. 



14 



REVTgX 3.0 



Cancio et al. Preprint, 2001 



Convergence studiesEBcH of a Slater- Jastrow trial wave- 
function for the Si atom reveal trends that we also ex- 
pect to be observed in the present study. The correlation 
hole averaged over the position of the reference electron 
should be quite close to that of the true ground state, 
with deeper minima and sharper maximum in proportion 
to the change in E c . The hole as a function of position 
in the system will in addition show subtle alterations in 
shape, while preserving basic qualitative features, for ex- 
ample the bond polarization in Fig. |]. 

The correlation hole can be expected to be more accu- 
rate at higher densities which carry greater weight in a 
variational optimization of the wavefunction; however the 
Slater- Jastrow wavefunction guarantees important condi- 
tions on the hole at any density, such as the cusp condi- 
tion and sum rules, so a reasonable estimate of the hole 
should be obtained even at quite low densities. Assuming 
that the generic trend at high density is to deepen the 
correlation energy of the hole by roughly 14%, we expect 
the overestimate of e c (r) by the DFT models, as shown 
for example in Fig. || to be reduced by between 10 and 
50% in the Si bond. Qualitative trends in this region 
will not be changed if the error in e c (r) does not vary 
dramatically with position. 

In general, the various errors in the VMC data (plane 
wave cutoffs, sum rule enforcement, and variational 
bias) are most significant for the angle-averaged n xc and 
exchange-correlaction energy per particle at high density. 
Even though estimated corrections are rather small, the 
agreement with the various DFT's studied is quite close 
and even small changes in the VMC data can be signifi- 
cant. Consequently, the errors in our data contribute to 
the general difficulty we find in assessing the high den- 
sity trends of the WDA and ADA. On the other hand, 
the difference between the DFT theories and the VMC 
are larger for the angle-averaged correlation hole. The 
qualitative trends in n c near bond center discussed in 
this paper are accordingly unaffected by error correc- 
tions though there are small differences in the quantita- 
tive value of e c (r) at the "bond-right" position. Finally, 
the improvement of ADA and WDA over LDA at low 
density is so dramatic that the observed error corrections 
to the VMC are small in comparison. 

VII. CONCLUSION 

We have carried out a detailed analysis of the 
exchange-correlation hole and energy per particle in the 
Si crystal, comparing various DFT models to accurate 
numerical data calculated with the VMC method. 

We find that the WDA and ADA help overcome the 
major defects of the local density approximation at low 
densities and especially in the pseudopotential atom core, 
where the rapid change in density relative to the length 
scale of the LDA hole dramatically affects the shape and 
range of the hole, and thus improve the fit with the VMC 



energy density. The remaining discrepancies are due to 
the failure of density averaging to provide significant in- 
formation at intermediate and high densities, where the 
inhomogeneity in the density is less severe, but the con- 
tribution to the energy density of subtle effects of the 
inhomogeneity is noticeable because of the higher den- 
sity. These discrepancies are, at least in principle, the 
result of the inflexibility of the scaling form of the pair 
correlation function used to fit the hole, which is insen- 
sitive to subtle variations in density. 

The detailed investigation of exchange and correlation 
at high density reveal the importance of orbital corre- 
lations in both cases. We find that the exchange hole 
is well described by a Wannier bonding state near the 
bond center, and that a "left-right" correlation or bond 
polarization plays an important role in the spherically- 
averaged correlation hole and has a noticeable effect on 
the correlation energy density. This bond polarization 
has the effect of countering the coarse-scaled deviation 
of the exchange hole from an isotropic form centered on 
the electron. As a result, the exchange-correlation hole 
and energy density is much more reliably fit by the DFT 
models studied than either exchange or correlation alone. 

Viewed in another sense, these orbital effects lead to se- 
rious defects in the correlation hole and energy of models 
derived from the homogeneous electron gas. The WDA 
and ADA conspicuously fail to improve on the LDA both 
in the system averaged energies E x and E c and in detail, 
particularly in the bonding region. Here, the nonlocal 
corrections to the LDA employed by these models fail to 
account for the physical features responsible for the most 
significant errors in the LDA model. 

Of all the methods considered in this paper, the GGA 
alone accurately predicts E x and E c for Si. Despite the 
lack of a direct comparison between our data and current 
GGA models for n xc ,c3 the assumptions made in con- 
structing the GGA exchange and correlation hole appear 
to be generally born out by our results, with a notable 
exception in the low density, large inhomogeneity limit 
for the correlation hole. Thus we expect that the accu- 
rate GGA results for E x and E c are due to an improved 
description of n x and n c on average, if not locally. Inter- 
estingly, the point by point trends of our data seem to 
be best characterized by the Laplacian of the density, in 
contrast to the approach taken by current gradient-only 
GGA models. Moreover, exchange and correlation taken 
alone show nonlocal qualitative features at high density 
that may prove difficult to describe in terms of any local 
density expansion. 

It would be of great value with respect to augmenting 
exact- exchange calculations to have a compact expression 
of the bond polarization features of the correlation hole 
in terms of localized atomic or bond-centered orbitals, 
possibly with a short-range correction in LDA or GGA. 
We have done preliminary calculations on modeling this 
bond polarization in terms of excitations of electrons into 
antibonding states. However a conveniently usable and 
compact form for the correlation hole at long range re- 
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mains to be found. 
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